#include <iostream>
#include "cblacs.h"
#include "commoninc.hpp"
#include "numvec.hpp"
#include "nummat.hpp"

using namespace std;
int main(int argc, char** argv){
  int myid, nprocs;
  int ctxt;
  int nprow, npcol;
  char roworder='R';
  char cA[] = "A", cZ[] = "Z";
  int myprow, mypcol;
  int MB, NB;
  int N;
  int lda, ldz;
  
  int DESCA[DLEN_], DESCZ[DLEN_];
  int IZERO = 0, IONE = 1, ISIX = 6;
  double DZERO = 0.0;
 
  int info; 
  DblNumMat      A, Z;
  DblNumVec      W;
  DblNumVec      work;
  IntNumVec      iwork;
  IntNumVec      ifail, iclustr;
  DblNumVec      gap;

  int lwork, liwork;

  char jobz     = 'V', 
       uplo     = 'L';

  /* MPI initialization */

  MPI_Init(&argc, &argv);
  MPI_Comm_rank(MPI_COMM_WORLD, &myid);
  MPI_Comm_size(MPI_COMM_WORLD, &nprocs);
  Cblacs_get(0,0, &ctxt);
  nprow = 1;
  npcol = 2;
  Cblacs_gridinit(&ctxt, &roworder, nprow, npcol);
  Cblacs_gridinfo(ctxt,  &nprow, &npcol, &myprow, &mypcol);
  
  /* Matrix initialization */
  MB = 2; NB = 2;
  N = 3;
  lda = 4; ldz = lda;
  
  descinit_(DESCA, &N, &N, &MB, &NB, &IZERO, &IZERO, &ctxt, &lda, &info);
  descinit_(DESCZ, &N, &N, &MB, &NB, &IZERO, &IZERO, &ctxt, &lda, &info);
  iA(info == 0);
  
  A.resize(lda, lda);
  setvalue(A, 0.0);
  if( myid == 0 ){
    A(0,0) = 1.0;
    A(1,0) = 1.0;
    A(2,0) = 1.0;
    A(0,1) = 1.0;
    A(1,1) = 2.0;
    A(2,1) = 3.0;
  }
  if( myid == 1 ){
    A(0,0) = 1.0;
    A(1,0) = 3.0;
    A(2,0) = 6.0;
  }

  work.resize(10*MB);
  pdlaprnt_(&N, &N, A.data(), &IONE, &IONE, DESCA, &IZERO, &IZERO, cA, &ISIX, &work[0]);

  Z.resize(ldz, ldz);
  W.resize(lda);
  setvalue(Z, 0.0);
  setvalue(W, 0.0);
  lwork = 1000*N;
  liwork = 1000*N;
  work.resize(lwork);
  iwork.resize(liwork);

  pdsyevd_(&jobz, &uplo, &N, A.data(), &IONE, &IONE, DESCA, 
	   W.data(), Z.data(), &IONE, &IONE, DESCZ, work.data(), &lwork, 
	   iwork.data(), &liwork, &info);
  iA(info == 0);


  if(myid == 0){
    cout << endl << W << endl;
  }
  work.resize(10*MB);
  pdlaprnt_(&N, &N, Z.data(), &IONE, &IONE, DESCZ, &IZERO, &IZERO, cZ, &ISIX, &work[0]);
  
  MPI_Finalize();
  return 0;
}
